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We propose a computational framework for the self-consistent dynamics of a microsphere system 
driven by a pulsed acoustic field in an ideal fluid. Our framework combines a molecular dynamics 
integrator describing the dynamics of the microsphere system with a time-dependent integral equation 
solver for the acoustic field that makes use of fields represented as surface expansions in spherical 
harmonic basis functions. The presented approach allows us to describe the inter-particle interaction 
induced by the field as well as the dynamics of trapping in counter-propagating acoustic pulses. 
The integral equation formulation leads to equations of motion for the microspheres describing the 
effect of non-dissipative drag forces. We show (1) that the field-induced interactions between the 
microspheres give rise to effective dipolar interactions, with effective dipoles defined by their velocities, 
and (2) that the dominant effect of an ultrasound pulse through a cloud of microspheres gives rise 
mainly to a translation of the system, though we also observe both expansion and contraction of the 
cloud determined by the initial system geometry. 

I. INTRODUCTION 

Computational approaches that employ an integral 
equation formalism to examine acoustic scattering from 
particles typically assume a static environment in which 
scatterers remain stationary. At present, a large body 
of work details such scattering problems [IHS]. While 
these stationary integral equation methods offer a large 
degree of accuracy in capturing the underlying physics, 
many problems of interest require a fully dynamical treat¬ 
ment. For instance, in biomedical physics, gas-filled micro¬ 
spheres exposed to ultrasonic beams have demonstrated 
effectiveness as a contrast imaging agent [1] and as drug 
delivery method [Sill], and Ding et al. have demonstrated 
their manipulation using acoustic tweezers in microflu¬ 
idic channels [7]. Moreover, composite materials consist¬ 
ing of colloidal in-fluid suspensions have peculiar sound 
propagation properties that can deviate from the ones of 
homogeneous liquids [8]. In each of these applications, 
the unconstrained motion of scatterers requires a self- 
consistent description of their dynamics in conjunction 
with a description of the acoustic field propagation. 

Here, we demonstrate the applicability of coupling par¬ 
ticle kinetics to a time-domain integral equation (TDIE) 
scattering framework to model rigid-sphere motion in¬ 
duced by a time-dependent acoustic potential. Specifi¬ 
cally, we consider the case of an acoustic pulse acting on 
microspheres that move in a fluid. Effective Langevin 
time-averaged radiation pressure forces EKH], which con- 
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sider the case of a steady radiation flux incident upon a 
body kept in static equilibrium, do not provide an appro¬ 
priate model in this case as they cannot accommodate 
inter-particle scattering effects. While many theoretical 
and computational descriptions of higher-order acoustic 
interactions exist [TT1 - H5] . few actually make use of com¬ 
puted fields to predict particle trajectories. As we consider 
only short-duration pulses, we refrain from time-averaging 
in favor of using a time-domain scattering formulation 
to explicitly calculate particle trajectories resulting from 
a prescribed pulse. By adopting a weakly-compressible 
potential formulation of the fluid media, our scalar wave 
problem inherits a number of similarities and solution tech¬ 
niques from scattering problems in electromagnetic the¬ 
ory, a topic previous works discuss extensively [niiiKn]. 
Moreover, our time-domain formulation readily allows the 
study of transient phenomena (such as acoustic tweezing); 
a convenience not shared with more common frequency 
domain approaches. 

We structure the remainder of this paper as follows: 
we first provide a formal mathematical description of 
the problem—including details on both the kinetic and 
field methods—followed by data obtained from various 
pulse and microsphere configurations, demonstrating both 
attractive and repulsive regimes suitable for subtle control 
of spherical systems in a homogeneous fluid. Einally, 
we offer concluding remarks on the effectiveness of the 
simulation as well as our thoughts on possible future 
extensions. 


II. CONTINUUM PROBLEM STATEMENT 

Consider a collection of N rigid, non-intersecting spher¬ 
ical scatterers (microspheres), each having radius a^, posi- 
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tion Xfe, and enclosing volume 14 C . The microspheres 
move in a homogeneous exterior fluid occupying Ve , where 
we denote the boundary of each microsphere as = i914 
and thus may ascribe to each an outward-pointing nor¬ 
mal rife(0,(/)), where 9 and 4> represent colatitude and 
azimuthal angles with respect to the local origin (micro¬ 
sphere center). We wish to investigate the reaction of 
the system to an incident acoustic pulse, thus the fluid 
carries a prescribed (band-limited) waveform through the 
microsphere system in which it interacts with each of 
the dflk according to the “sound-hard” regime presented 
in m- The incident acoustic pulse, in combination with 
the acoustic field scattered from each microsphere and 
the hydrodynamic field induced by the relative velocity of 
each microsphere, acts as a perturbation to the initially 
at-rest uniform ideal fluid [HUS]. We consider here the 
linear regime, in which the perturbation induced by the 
acoustic and aerodynamic contribution remain sufficiently 
small so that the velocity field v(x,t) satishes the condi¬ 
tion |v(x, t)| ^ Cs, where Cg represents the speed of sound 
in the fluid. In this limit, the velocity potential, defined 


(/7(x,t) = v?i„c(x,t) -f ^ f dt' f dA 

i—Q V 


by v(x, t) = V(^(x, t), satisfies the scalar wave equation: 

and we may express the pressure perturbation at any 
point in the exterior medium as 

dip{x,t) 

p{x,t) = -po ——, (2) 

where po denotes the equilibrium density of the fluid. 
Rigidity of the fife necessarily prescribes boundary condi¬ 
tions on the normal velocity components at each interface, 
namely. 


dpjyiA) 

drik 


dxfe 





(3) 


where Xfe represents the center-of-mass coordinate for the 
fc**' microsphere. 

Using these relations, we apply the Kirchoff-Helmholtz 
theorem to define the following system of integral equa¬ 
tions, 


9Gr.(x, t; x',t') 
^fife 


Gr(x, t. 


x',t') 


Sftfe / ’ 


(4) 


where Gr(x,t; x',t') denotes the Green’s function for a 
retarded potential. 


Gr{x,t; x',t') 


6{t - f' - |x - x'l/cg) 
47r|x — x'l 


(5) 


If the system remains localized to a region with small 
dimensions when compared to the wavelength of sound, 
retardation effects become negligible and we may instead 
use the Laplace-kernel Green’s function. 


G(x,x') 


1 

47r|x — x'l 


( 6 ) 


To ease notation, we define the following two integral 
operators, 


5fe[(/5(x e flfe(t),t)] = f d4G(x,x') (7a) 

"'Ofc(t) 

X>fe[(/3(x e flfe(t),t)] = f dAip{x',t) 5fl^G(x,x'), (7b) 

jQk{t) 

reducing Eq. @ to 

7V-1 

(p{x, t) = (^inc + y] [Vk - Sk^ [(^(x e flfe(t), <)]. (8) 

In solving Eq. Q, we obtain the velocity potential ev¬ 
erywhere for a given time without retarded scattered fields. 


For the incident pulse, ^pi■ac^ we consider superpositions 
of wave packets of the form 

VJinc(x,t) =PoCOs(wot-k-x)e“(“'’‘“^’") (9) 

Finally, the variation in pressure (and thus p) over each 
of the fife necessarily propels each microsphere according 
to the equation of motion 


ruk 


d^xfc 

dt^ 


Pq 


n^{t) dt 


( 10 ) 


III. DISCRETIZATION OF THE INTEGRAL 
EQUATIONS 

To solve the integral equation scattering problem, we 
begin by discretizing our field in both space and time. As 
we have restricted our particles to completely spherical 
geometries, the spherical harmonics, dehned by 


YiM) = ( 11 ) 

give simple eigenfunctions of the operators in Eq. (7b). 
As a result, they lend themselves well to an expansion of 
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FIG. 1. Coordinate notation. 


Lp on the surface of each microsphere with respect to the 
microsphere’s center, 


t>0 \m\<i 


( 12 ) 


By considering Eq. ([^ and expressing the local velocity 
potential at each of the 12^ as a linear combination of 
spherical harmonics, we have a complete representation 
of the body force acting on each microsphere. 


FLdyW = - [ dSp(x e 
J Q,k(t) 


cfi(t) - Cti(i) 


+ c'l-My- (13) 


— Po\ 


due to the orthogonality of dipole terms with the rest of 
the multipoles. 

The problem then becomes one of solving a system of 
linear equations that we may compactly represent as 




(14) 


with the overbar denoting a matrix quantity. We define 
the elements of iF as projections of the incident field onto 
local spherical harmonics, 

FL = IdA Yl^{e, 0)(^inc(x, t), (15) 

jQk{t) 

and detail for two cases: j = k and j ^ k. In the 


instances where j = k, Eq. (7b I propagates effects of the 


interaction through to every point on a surface sharing a 
coordinate system with the original, thus the harmonics 
remain orthogonal and 


■P'33 


1 


2^ + 1 


(16) 


after exploiting the well-known expansion theorem for 
Eq. ([^, 

G(x,x') = ^ ^^^y,„(0,(/.)F;„(0',</)') (17) 

where r< = min(|x|, |x'|) and r> = max(|x|, |x'|). A 
description of the off-diagonal terms where j ^ k proceeds 
much the same way, though the surface expansions no 
longer share a local origin, complicating the projections. 
Translation operators for the spherical harmonics EHl 
[2T] allow analytic expressions for these matrix elements, 
though we eschew such operators in favor of numerical 
integration for speed. 

Thus, at every timestep of the simulation, the algo¬ 
rithm proceeds as follows: (i) project the incident pulse 
and surface velocities onto local expansions of spherical 
harmonics, (ii) propagate scattering effects through space 
by inverting the operators in Eq. (§, (hi) project these 
scattered fields onto local spherical harmonics to give a to¬ 
tal representation of tp on each surface, and (iv) move each 
microsphere according to Eq. (10) & advance t ^ t + At. 
For rigid microspheres only i = 1 terms contribute to 
center-of-mass motion, thus we use only the Cim coeffi¬ 
cients in evolving Eq. (10). 

The inversion in step (ii) above requires some care; by 
simply inverting the entire propagation operator, V — S, 
to give a single surface pressure, Eq. (10) reduces to a 
differential equation of the form 


Xfc = /(t,Xfc,Xfc). 


(18) 


This presents a number of irregularities with conventional 
integration schemes and will rapidly diverge towards ±oo 
due to the additional x^ on the right if implemented 
naively. To remedy this, we note that S serves to produce 
only a reaction or drag term on each microsphere that im¬ 
pedes motion. By maintaining quantities for the inversion 
of T) and S separately, we remove the explicit dependence 
on Xj; by introducing a linear coefhcient in the form of an 
additional mass term—given by the x^-dependent con¬ 
tribution in the single-layer S operator—when solving 
Eq. (pi). 


IV. ANALYTIC RESULTS 
A. Single microsphere solution 

As an example, consider a single sphere of density ps 
and radius a. Taking ka 1, we may approximate Eq. (|^ 
as (finci'^A) = yo{t)z and we wish to find the response 
velocity of the sphere, u, in terms of the field velocity 
V = V:^inc- It follows that the expansion of ^^inc contains 
only ^ = \ terms, thus 

(/Sine = iio(t)acos(d) 


^mm' 


(19) 
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on the surface of the sphere. Similarly, from Eq. (|^, 
Sfltp = u • n 

= Uzacos{0) ( 20 ) 

due to the symmetries present in x and y. As a result, 


f ~ J dS''v5(x')9fi/G'(x, x') = 

voacos{9) — J dS'auz cos(0) G(x, x'). 


( 21 ) 


and it becomes apparent that only ^ = 1, m = 0 terms in 
Eq. @ remain after integrating. Consequently, the field 
becomes 


V3(x,t) = Uo(t)|x| + 


a^{vo{t) - Uz) 


2|x|' 


cos(6») (22) 


outside the microsphere and 


Vj(x e fl,t) = ( ^uo(t) - ]^Uz )acos(6») (23) 


Ul U2 



FIG. 2. Perpendicular configuration. 


B. Low-order interactions 

We now consider two identical microspheres arranged 
perpendicularly to an incident waveform as in Fig. 
Within the Born approximation, we may take Eq. as 
the incident field and use it in place of the total field on 
the right-hand side of Eq. Q, assuming negligible contri¬ 
butions from scattering. In doing so, the field everywhere 
becomes 


on its surface. From this we conclude the total velocity 
potential in the fluid arises from a surface-scattering term 
alongside a term describing the transfer of momentum 
from the moving microsphere to the fluid. 

Using Eq. (10), we may then write the equation of 
motion for the system as 


PsViiz 


PoV 



(24) 


where V = 47ra^/3 gives the volume of the microsphere. 
The transfer of momentum from the moving microsphere 
to the fluid becomes a reaction force of the fluid due to the 
sphere. Landau & Lifshitz |19j initially derived this non- 
dissipative drag force by way of momentum and energy 
conservation. Note that this drag force presents only in 
the case of accelerated motion of the microsphere and we 
may recast its effect in the form of a virtual mass that 
includes a contribution due to the mass of the displaced 
fluid. 



This expression leads to a simple relation linking Uz(t) 
and vo(t) provided the velocity does not remain constant 
and that the sphere does not move in the absence of the 
held: 


_ 3po 

uo Po + 2ps 


(26) 


The idea of a virtual mass for the accelerated motion 
of a single sphere in an ideal fluid readily generalizes to 
the case of a moving collection of mutually-interacting 
spheres. Through this, we may compute the dynamics 
of each microsphere in the group, taking into account 
the effect of the momentum exchange between the fluid 
and the microspheres, resulting in both drag and inter¬ 
particle forces in addition to the displacement caused by 
the driving acoustic field. 


, , , , a^ cos(0i) r T 

Vj(x,t) =vo{t)z+ -;—;:T 2 UO-wiJ 

4 |x- di 2 / 2 | 


cos( 02 ) 


3 |x-t-di 2 / 2 |' 


• [vo - U 2 ] . 


(27) 


By inserting this into Eq. (101 for xi, we have 

miUi • z = 27rpoa^ J cos^ 0ia^ ~ d(cos 61 ) + 


Po 


a° Vo — U 2 
\ 3 |x-di2| 


■2 cos 01 cos 02 d^i d(cos 0i). 


Writing 


cos 02 = 


cos 01 


^12 


y(l-^sin0i)' + (^cos0i)' 


(28) 

(29) 


and noting Ui = U 2 = due to symmetry in the initial 
configuration, we may expand Eq. (28) in a/di 2 to give 


PsUs = Pol ) -f 


po{vo - Us) f a 


In the limit of di 2 —> 00 , this becomes 

u^ _ 4po 
1^0 Po + 3ps 




( 31 ) 


By considering negligible scattered fields at the surface 
of each microsphere, we qualitatively recover Eq. (27) 
with different coefficients arising only from the Born ap¬ 
proximation. Moreover, the additional interaction term 
in Eq. (31) scales as |d.y | a behavior anticipated from 
the dipolar nature of Eq. (22). 
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Quantity 

Symbol 

Value 

Sound speed 

Cs 

1500ms“^ 

Microsphere radius 


1 pm 

Density (exterior) 

po 

1000 kgm“^ 

Density (interior) 

ps 

Ikgm"® 

Pulse amplitude 

Po 

0.05 m^s”^ 

Center frequency 

fo 

0.5 MHz to 20 MHz 

Pulse duration (st. dev.) 

a 

7 ps to 24 ps 


TABLE I. Typical simulation parameters. 



FIG. 3. Translation of a single microsphere interacting with 
an incident pulse (/o = 0.5MHz, a = 7ps). Microsplieres 
interacting with the pulse translate a finite distance along k 
due to the Gaussian envelope in Eq. 


V. NUMERICAL RESULTS 


Here we present a series of numerically-solved systems 
to illustrate the utility of the method in investigating 
acoustic phenomena. We perform simulations of one- 
and two-particle/pulse systems to determine the principal 
particle-field and particle-particle interactions, followed 
by simulations of larger assemblages of spheres to inves¬ 
tigate group phenomena and effects in systems without 
symmetry. Unless otherwise stated. Table |T] gives the sim¬ 
ulation parameters for each of the following simulations; 
as our interests lie in hydrodynamic applications, we use 
material parameters characteristic of water to define our 
external fluid medium. Similarly, we consider here the 
case of gas-filled microspheres [1], and therefore set their 
density much smaller than that of the exterior medium. 
The acoustic pulses lie in the ultrasonic regime, and the 
chosen frequency of 20 MHz corresponds to that of typical 
applications in acoustic microscopy. 


Pressure gradient, z ■ Vp (relative scal^ 


100 


-50 


-100 



60 80 100 120 140 

Time (ps) 


FIG. 4. (Color online) Confinement of non-interacting spheres 
to planes; identical counter-propagating pulses (/o = 20 MHz, 
a = 23.8 ps) initially displaced along z tend to align objects in 
VP = 0 planes at A/2 intervals. Field & trajectories sampled 
every 30 timesteps and smoothed with a 16-sample windowed 
average. 


A. Single microspheres 


Figure gives the trajectory of a single microsphere 
initially at rest under the effects of an incident Gaussian 
pulse. Under the linear and ideal fluid approximations 
and absent the Gaussian envelope in Eq. § , the micro- 
sphere merely oscillates about its origin in accordance 
with Eq. (26). In the pulsed case, however, the variation 
in pressure imposed by the finite value of k modifies the 
system dynamics to yield a net translation of each micro- 
sphere. Note that the regime considered here produces 
no net transfer of momentum between the acoustic field 
and the microsphere—a consequence of the ideal fluid. 


Figure depicts smoothed results of 128 trajectories 
corresponding to single microspheres initially spaced along 
z and excited by identical counter-propagating pulses. By 
taking the width of each pulse much greater than the 
radius of each microsphere, the two pulses reproduce the 
effects of interfering standing waves. The confinement 
occurs at VP = 0 (nodal) planes where the net force 
on each microsphere vanishes. The half-wavelength as¬ 
sociated with the dominant pulse frequency gives the 
separation between neighboring planes. 


Finally, Fig. shows the relative velocity potential 
near a single microsphere; given a surface expansion of 
we may compute the potential everywhere through 
application of Eq. (§. As predicted by Eq. ( |^ , this 
field greatly resembles that of a pointlike “velocity dipole” 
with Vg acting as a dipole moment. 

The simulations described thusfar demonstrate precise 
acoustic control; through careful application of the inci- 
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FIG. 5. (Color online) Calculated isopotential contours 
near a lone microsphere. Red and blue colorations represent 
regions of positive and negative potential. The motion of each 
microsphere through the background medium serves primarily 
to produce a dipolar field of velocity potential with Vs serving 
as the sphere’s dipole moment. 


dent field parameters, we may induce a (finite, given a 
finite pulse) translation along the principal k-vector with 
a large degree of accuracy in the overall displacement. 
In addition, the application of multiple pulses serves to 
confine microspheres to highly localized regions in space, 
offering a self-consistent model of acoustic tweezing. 


B. Many-particle simulations 


We now turn our attention to collections of mutually 
interacting microspheres. To quantify the effects of scat¬ 
tering, we first decouple scattering forces from the incident 
pulse by arranging two microspheres perpendicularly to 
the pulse’s k-vector. Figure gives results for such a 
simulation where we plot the relative change in velocity 
as compared with the single-particle simulation. 


^|Vmax| — Cnax(|vuouble(0 ^single (^) |) ■ 


In principle, describing quantities found from a complete 
simulation as a function of initial separation could obfus¬ 
cate scaling data considerably; forces arising from scatter¬ 
ing could alter the geometry of the system. In practice, 
however, the perpendicular configuration used here gives 
scattering forces that only influence the motion along k. 
Consequently, Av oc z and the microspheres’ initial sepa¬ 
ration remains a good estimator of scaling behavior. We 
see in Fig. that the radii data scale as and the separa¬ 
tion data exhibit strong |di 2 | ^ scaling, again indicating 
a dominant dipolar interaction between microspheres as 
shown by Ilinskii et al. in 2007 m and predicted by 
Eq. (|22l). 


FIG. 6. (Color online) Scaling behavior of two microspheres 
arranged perpendicularly to an incident pulse for various 
radii and initial separations. The (a, •) symbols on each 
axis denote data associated with that axis. The a follow a 
regression of A|v|j = 0.250 754d)-2^ “““^^, and the • follow 
A|vinax|r = 3.133 28 X These trends strongly 

indicate dominant dipolar interactions between microspheres. 



FIG. 7. (Color online) Isosurfaces of velocity potential 
(arb. units) calculated by evaluating the 5 and ID terms in 
Eq. (|^ for a A = 16 particle simulation. Red, blue, and yellow 
surfaces denote regions of positive, negative, and zero potential, 
with holes appearing due to intersections with the bounding 
box. The inset box shows the three dimensional arrangement 
of the microspheres superimposed with their velocity vectors, 
as well as several positive and negative potential isosurfaces. 
Rendered with Visit [22]. 
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FIG. 8. (Color online) Fractional change in the volume of 20 
randomly-initialized microsphere clouds subject to the same 
incident pulse, smoothed with a 128-sample moving average. 
Positive and negative values denote expansion and contraction. 
a = 1.5 cm. 


Finally, we consider the dynamics of large {N = 16) 
clouds of microspheres. For each simulation, we generate 
a collection of microspheres initialized with zero veloc¬ 
ity and random positions within a 10 pm ball subject to 
a minimum-separation constraint to prevent collisions. 
Figure shows a snapshot of the velocity potential iso¬ 
surfaces calculated in one such simulation. Even with 
mutual interactions, the shape of each isosurface remains 
consistent with the presence of a dipolar field oriented 
along the microspheres’ velocity. Again, due to the local¬ 
ization assumption used to justify Eq. (§, each system 
predominantly translates a finite distance in accordance 
with the results found for a single microsphere in Fig. 

To quantify small changes in the geometry of a system, 
we compute 14,, the volume of the convex hull containing 
each microsphere, at every timestep in the simulation [2dj . 
Figure]^ shows the fractional change in the hull volume. 


^Vh 


Vh{t)-Vh{0) 

VhiO) 


(33) 


for 20 such systems after smoothing with a weighted 
moving average. Curves ending above and below zero 
indicate larger and smaller hull volumes (system expansion 


and contraction). We note from Fig. |^a greater tendency 
for random clouds to expand; the effective dipole-dipole 
interaction between particles with T k gives purely 
repulsive forces, while the interaction between particles 
with dij II k gives both repulsive and attractive effects 
depending on a and the relative phase of the oscillating 
microsphere velocities. 

VI. CONCLUSIONS 


This work lends a novel, fine-grained approach to the 
study of acoustic response via integral equation meth¬ 
ods. By considering a potential representation in terms 
of spherical harmonics on the surfaces of microspheres 
coupled to a standard molecular dynamics scheme, we 
obtain a description of the microspheres’ dynamics under 
the effect of ultrasound pulses without resorting to time- 
average approximations, though the confined microsphere 
geometries under consideration allow us to neglect small 
effects arising from time-delays in scattering. We have 
shown that the net effect of an ultrasound pulse on a sin¬ 
gle microsphere consists of a translation that we can tune 
through careful control of pulse parameters. Additionally, 
systems with multiple incident waveforms tend to confine 
microspheres to nodes in the pressure field governed by 
acoustic interference. Finally, in the dynamics of systems 
with many microspheres, we have observed the effect of 
weak inter-particle transient effects induced by the driving 
acoustic pulse. These effects can produce both expansion 
and contraction of a cloud of microspheres, in addition 
to the overall translation. 

Prior work in this area [24] [2^ makes use of deformable 
bubble boundaries about fixed locations. Incorporation 
of these methodologies to our theoretical model natu¬ 
rally offers possibilities for future research, as does the 
addition of retardation effects. Additionally, we expect a 
straightforward approach to experimental confirmation 
of the results presented here. Optical tracking of tracer 
particles has demonstrated its effectiveness in similar 
fluid-trajectory studies and would readily adapt to track 
physical analogues of our theoretical microspheres. 
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